MAE0552 - Introdução à Teoria da Informação
Projeto final: Um estudo sobre o sequenciamento genético e sua relação com características dos indivíduos
| Ana Luisa Pinheiro | 11810407 | |
| Ayrton Amaral | 11288131 | |
| Bruno Groper Morbin | 11809875 | |
| Caio Febronio | 11811482 |
Instituto de Matemática e Estatística - Universidade de São Paulo | Julho, 2023
Para cada indivíduo pertencente ao estudo, foram amostrados genes individuais e realizada a medição da intensidade da luz em cima de cada um deles. O banco de dados disponibilizará os valores dessa intensidade individualmente para cada gene de cada indivíduo.
Com base no banco de dados, queremos verificar se existe alguma relação entre esses valores da intensidade da luz em cima de cada gene com características pessoais de cada indivíduo.
Para a resolução do problema, agruparemos os genes pelo valor da intensidade de luz e utilizaremos a informação mútua vista em aula para checarmos se os genes introduzidos em cada cluster possuem distribuição parecida ou não. Após essa checagem, será feita uma distribuição única da intensidade de luz para cada cluster e a mesma será utilizada para caracterizar cada conjunto de cluster em relação às características do indivíduo.
# Carregando pacotes
library(tidyverse)
library(dplyr)
library(cluster)
library(infotheo)
library(stringr)
# library(ggpubr)
library(patchwork)load("glioma.RData") # geneInfo ; gliomaGSE52009 ; targetInfoGlioma
glioma <- gliomaGSE52009[1:50,]; as.data.frame(glioma)info <- targetInfoGlioma; rownames(info) <- NULL; info |> select(colnames(info[,-1]),FileName)str(info[,-1]) # ignorando a coluna FileName'data.frame': 120 obs. of 7 variables:
$ gender : chr "female" "female" "female" "unknown" ...
$ age : num 37 61 37 -100 31 42 42 38 32 40 ...
$ geoAccession: chr "GSM1257398" "GSM1257399" "GSM1257400" "GSM1257401" ...
$ sampleInfo : chr "Astrocytoma.a-0253" "Astrocytoma.a-0258" "Astrocytoma.a-0281" "Astrocytoma.a-0285" ...
$ diagnostic : chr "astrocytoma" "astrocytoma" "astrocytoma" "astrocytoma" ...
$ datasetId : chr "Glioma52009" "Glioma52009" "Glioma52009" "Glioma52009" ...
$ tissue : chr "brain" "brain" "brain" "brain" ...
info$gender <- factor(info$gender); levels(info$gender) "female" "male" "unknown"
ggplot(info, aes(x=reorder(gender, -table(gender)[gender])))+
geom_bar(aes(fill=gender), color="transparent")+
scale_fill_manual(values=c(male="#3E67A3",female="#A34336",unknown="#7F8F85"))+
geom_text(stat = 'count', aes(label = paste0(round((after_stat(count)/sum(after_stat(count)))*100), "%")), vjust = -0.5) +
scale_y_continuous(expand = expansion(mult=c(0,.2)))+
labs(x=NULL,y=NULL, title="Gênero")+
guides(fill="none")info$diagnostic <- factor(info$diagnostic); levels(info$diagnostic) "anaplastic.astrocytoma" "anaplastic.oligodendrocytoma"
"anaplastic.oligodendroglioma" "astrocytoma"
"glioblastoma" "oligodendroglioma"
ggplot(info, aes(x=reorder(diagnostic, -table(diagnostic)[diagnostic])))+
geom_bar(aes(fill=diagnostic), color="transparent")+
scale_fill_grey(end = 0.9, start=.5)+
geom_text(stat = 'count', aes(label = paste0(round((after_stat(count)/sum(after_stat(count)))*100), "%")), vjust = -0.5) +
scale_y_continuous(expand = expansion(mult=c(0,.15)))+
labs(x=NULL,y=NULL, title="Diagnóstico")+
theme(axis.text.x = element_text(angle = 45,hjust = 1, vjust=1))+
guides(fill="none")range(info$age) -100 70
cat(paste0(sum(info$age<=0)," entradas inválidas para idade => ",round(sum(info$age<=0)/nrow(info)*100,2), "% da amostra")) # Porcentagem de dados inválidos para idade17 entradas inválidas para idade => 14.17% da amostra
range(info[which(info$age>0),]$age) 17 70
ggplot(subset(info[which(info$age>0),]), aes(x = age, y = after_stat(density))) +
geom_histogram(aes(y = ..density..), fill = "skyblue", color = "#0c0c0c", binwidth = 5, alpha = 0.9) +
geom_density(color = "cyan", linetype = "solid", linewidth = 1, fill="transparent") +
labs(x = "Idade", y = "Frequência relativa", subtitle ="(somente dados válidos)", title="Distribuição de idade na amostra")+
scale_y_continuous(expand = expansion(mult=c(0,.20)))+
scale_x_continuous(n.breaks = 20)unique(info$datasetId) "Glioma52009"
unique(info$tissue) "brain"
Aplica-se um agrupamento hierárquico para simplificar conjunto de genes que apresentam mesma intensidade de luminosidade para cada indivíduo. Dessa forma, e levando em consideração que os dados são coletados sob a mesma medida, opta-se por não padronizá-los. Caso fossem padronizados, algumas intensidades não expressivas poderiam ser conectadas com outras mais expressivas originalmente (quando considerados os valores não padronizados).
custom_dist <- function(x, y) exp((x-y)**2) # penaliza mais para grandes erros
d <- proxy::dist(glioma[,1], method=custom_dist)
for(i in 2:ncol(glioma)){
d <- d + proxy::dist(glioma[,i], method=custom_dist)
}Para a junção dos clusters, tem-se algumas ordens de hierarquização, tais quais:
Single linkage: É a distância mais curta entre quaisquer dois pontos nos dois grupos.
Complete linkage: É o oposto do single linkage. É a distância mais longa entre quaisquer dois pontos nos dois grupos.
Average linkage: É a média das distâncias entre cada ponto de um grupo para todos os pontos do outro grupo.
Centroid linkage: A distância entre o ponto central de um grupo ao ponto central do outro grupo.
Ward’s linkage: Uma combinação dos métodos average e centroid. A variância dentro do grupo é calculada determinando o ponto central do grupo e a distância das observações em relação ao centro. Ao tentar mesclar dois grupos, a variância é calculada entre os grupos e os grupos são mesclados se sua variância for menor em comparação com outras combinações.
hc<-hclust(d, method = "ward.D") # Hierarchical Cluster com linkage por Ward.D
num_cluster_init <- 20; init_clusters = cutree(hc, num_cluster_init) # pegando os clusters iniciais # Exemplo do um cluster selecionado inicialmente
as.data.frame(glioma[which(init_clusters==1),])genes.clusters <- list()
for(clust in 1:num_cluster_init){
gp <- glioma[which(init_clusters==clust),] # selecionando o cluster
if(is.null(nrow(gp))) {
next
} else{
bounds <- apply(gp,MARGIN = 2, FUN = function(x) quantile(x,c(0.15,0.85))) # definindo os intervalos interquantil para cada indivíduo nesse conjunto de genes
# identificando cada intensidade dentro do cluster se está dentro do interquantil do indíviduo
in_bound <- as.data.frame(lapply(1:ncol(gp), FUN= function(i) {as.numeric(between(gp[,i],bounds[1,i],bounds[2,i]))})); colnames(in_bound) <- colnames(bounds); rownames(in_bound) <- rownames(gp)
# manter no cluster apenas os genes que se apresentaram bastante no todo da amostra dentro do intenquartil da intensidade de cada indivíduo
select_genes <- apply(in_bound, MARGIN = 1, function(x) sum(x)) >= 0.7*ncol(glioma)
select_names <- names(which(select_genes==T))
if(length(select_names)>1){
genes.clusters <- append(genes.clusters,list(select_names))
}
}
}
for(i in 1:length(genes.clusters)){
if(i==1) cat("--- Genes:\n")
cat(paste0("\nCluster ",i,":\n"))
cat(genes.clusters[[i]],sep=" | ")
cat("\n")
}--- Genes:
Cluster 1:
10 | 10000 | 10001 | 10009 | 100128782
Cluster 2:
10002 | 10003 | 10004 | 100124700 | 100128236 | 100128327 | 100128979 | 100129033 | 100129066 | 100129637
Cluster 3:
100127886 | 100128843 | 100129455 | 100129726
Agora que os genes já foram agrupados por similaridade, visualiza-se o valor da informação mútua entre e intra os grupos com os genes selecionados. Para isso, será necessário discretizar as entradas para então obter as distribuições empirícas de cada gene.
https://search.r-project.org/CRAN/refmans/infotheo/html/multiinformation.html
breaks <- seq(floor(min(glioma)), ceiling(max(glioma)), by=1) # definindo intervalos para discretizar
mutualinfo_clust <- list()
for(clust in 1:length(genes.clusters)){
cluster_disc_temp <- apply(glioma[genes.clusters[[clust]],], MARGIN = c(1,2),function(x) cut(x, breaks)) # discretizando cluster
mutualinfo_clust <- append(mutualinfo_clust,multiinformation(t(cluster_disc_temp))) # Correlação Total (em nats - quando usado logaritmo natural para entropia)
}
for(i in 1:length(genes.clusters)){
if(i==1) cat("--- Informação Mútua (em nats):\n")
cat(paste0("\nCluster ",i,":\t ",format(round(mutualinfo_clust[[i]],2),width = 6, nsmall = 2)))
}--- Informação Mútua (em nats):
Cluster 1: 0.56
Cluster 2: 0.22
Cluster 3: 0.98
ind_bound <- max(quantile(unlist(mutualinfo_clust),0.25)[[1]],1); ind_bound # métrica para clusters consistentes 1
clusters_consistentes <- which(mutualinfo_clust<ind_bound); clusters_consistentes # Clusters consistentes (com distribuição de intensidade luminosa similar) 1 2 3
df_consistentes <- data.frame()
for(clust in clusters_consistentes){
# definindo a mediana das intensidades luminosas para cada indivíduo da amostra considerando todos genes do cluster
cluster_dist <- apply(glioma[genes.clusters[[clust]],], MARGIN=2, FUN= function(x) median(x))
df_consistentes <- rbind(df_consistentes, data.frame(t(cluster_dist), row.names = clust))
}
# distribuições das intensidades lumninosas unificada por cluster
df_consistentesTem-se 3 características de cada indívíduo disponibilizadas conforme visto na seção de apresentação do conjunto de dados: gênero, idade e diagnóstico. Com isso, verifica-se se há distinção da distribuição de cada nível de tais características por cluster de genes consistentes obtidos anteriormente.
Contudo, inicialmente, é realizado um breve tratamento nos dados: para a análise dos níveis de gênero, desconsidera-se as entradas como “unknow”; além disso, categoriza-se a idade, e despreza-se os dados com entrada inválida como idade negativa; por fim, para a informação sobre diagnóstico, separa-se em 6 diagnósticos já apresentados, e em dois tipos - anaplástico e plástico.
Dessa forma, para cada cluster, tem-se a visualização da distribuição da luminosidade para esses diferentes grupos, e espera-se identificar clusters que discriminizam bem alguma característica.
valid_gender <- info |> filter(gender != "unknown") |> select(id = geoAccession, gender)
valid_age <- info |> filter(age > 0) |> select(id = geoAccession, age)
valid_age$cat_age <- cut(valid_age$age, breaks = c(0,7,14,20,30,40,50,60,70,80,100))
valid_diag <- info |> select(id = geoAccession, diagnostic) |>
mutate(diag_type = case_when(
str_detect(diagnostic, "anaplastic") ~ "anaplastic",
TRUE ~ "plastic"
)
)for(clust in clusters_consistentes){
meio_p1 <- sum(range(df_consistentes[as.character(clust),valid_gender$id]))/2
med_p1 <- median(unlist(df_consistentes[as.character(clust),valid_gender$id]))
leg_dir_p1 <- ifelse(med_p1<meio_p1, 1, 0)
p1 <- data.frame(int_luz = unlist(df_consistentes[as.character(clust),valid_gender$id]), gender = valid_gender$gender) |>
ggplot(aes(x = int_luz, fill = factor(gender), color=factor(gender))) +
geom_density(alpha = 0.2) +
scale_fill_manual(values = c(male="#1775ff",female="#ff371c"), name=NULL)+
scale_color_manual(values = c(male="#1775ff",female="#ff371c"), name=NULL)+
scale_y_continuous(expand = expansion(mult=c(0,.10)))+
scale_x_continuous(expand= c(0,0))+
theme(legend.position = c(leg_dir_p1,1), legend.justification = c(leg_dir_p1,1), legend.background = element_rect(fill = "#00000050"), plot.title = element_text(hjust=0), legend.key = element_rect(color="transparent", fill="transparent"))+
labs(x= "Intensidade luminosa", y="Densidade", title="Gênero")
meio_p2 <- sum(range(df_consistentes[as.character(clust),valid_age$id]))/2
med_p2 <- median(unlist(df_consistentes[as.character(clust),valid_age$id]))
leg_dir_p2 <- ifelse(med_p2<meio_p2, 1, 0)
p2 <- data.frame(int_luz = unlist(df_consistentes[as.character(clust),valid_age$id]), age = valid_age$cat_age) |>
ggplot(aes(x = int_luz, fill = factor(age), color=factor(age))) +
geom_density(alpha = 0.2)+
scale_fill_discrete(name=NULL)+
scale_color_discrete(name=NULL)+
scale_y_continuous(expand = expansion(mult=c(0,.10)))+
scale_x_continuous(expand= c(0,0))+
theme(legend.position = c(leg_dir_p2,1), legend.justification = c(leg_dir_p2,1), legend.background = element_rect(fill = "#00000050"), plot.title = element_text(hjust=0), legend.key = element_rect(color="transparent", fill="transparent"))+
labs(x= "Intensidade luminosa", y="Densidade", title="Faixa etária")
meio_p34 <- sum(range(df_consistentes[as.character(clust),valid_diag$id]))/2
med_p34 <- median(unlist(df_consistentes[as.character(clust),valid_diag$id]))
leg_dir_p34 <- ifelse(med_p34<meio_p34, 1, 0)
p3 <- data.frame(int_luz = unlist(df_consistentes[as.character(clust),valid_diag$id]), diagnostic = valid_diag$diagnostic) |>
ggplot(aes(x = int_luz, fill = factor(diagnostic), color=factor(diagnostic))) +
geom_density(alpha = 0.2) +
guides(fill=guide_legend(title=NULL), color=guide_legend(title=NULL))+
scale_y_continuous(expand = expansion(mult=c(0,.10)))+
scale_x_continuous(expand= c(0,0))+
theme(legend.position = c(leg_dir_p34,1), legend.justification = c(leg_dir_p34,1), legend.background = element_rect(fill = "#00000050"), plot.title = element_text(hjust=0), legend.key = element_rect(color="transparent", fill="transparent"))+
labs(x= "Intensidade luminosa", y="Densidade", title="Diagnóstico")
p4 <- data.frame(int_luz = unlist(df_consistentes[as.character(clust),valid_diag$id]), diagnostic = valid_diag$diag_type) |>
ggplot(aes(x = int_luz, fill = factor(diagnostic), color=factor(diagnostic))) +
geom_density(alpha = 0.2) +
scale_fill_manual(values=c("anaplastic"="#5ed19f","plastic"="#c753f5"))+
scale_color_manual(values=c("anaplastic"="#5ed19f","plastic"="#c753f5"))+
guides(fill=guide_legend(title=NULL), color=guide_legend(title=NULL))+
scale_y_continuous(expand = expansion(mult=c(0,.10)))+
scale_x_continuous(expand= c(0,0))+
theme(legend.position = c(leg_dir_p34,1), legend.justification = c(leg_dir_p34,1), legend.background = element_rect(fill = "#00000050"), plot.title = element_text(hjust=0), legend.key = element_rect(color="transparent", fill="transparent"))+
labs(x= "Intensidade luminosa", y="Densidade", title="Tipo de diagnóstico")
ptot1 <- p1+p2+plot_layout(ncol = 2) & theme(
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.background = element_rect(fill="#0c0c0c", color="#0c0c0c"),
text = element_text(size=8)
)
ptot2 <- p3+p4+plot_layout(ncol = 2) & theme(
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.background = element_rect(fill="#0c0c0c", color="#0c0c0c"),
text = element_text(size=8)
)
cat("<pre>--- Cluster ", clust, ":</pre>","\n")
cat("\n")
print(ptot1)
print(ptot2)
cat("\n\n")
}--- Cluster 1 :
--- Cluster 2 :
--- Cluster 3 :